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I. INTRODUCTION 



A quantitative prediction of fluid flow, sound propagation, or chemical transport in strongly correlated disordered 
media, such as sedimentary rock, frequently employs representative microscopic models of the microstructure as input. 
A large number of microscopic models have been proposed over the years to represent the microstructure of porous 
media 

Microscopic models do not reproduce the exact microstructure of the medium at hand, but are based on the idea that 
the experimental sample is a representative realization drawn from a statistical ensemble of similar microstructures. 
Hence it is necessary to have methods for distinguishing microstructures quantitatively [jl6|-|l9|. This is particularly 
important for attempts to generate porous microstructures in an automatic computerized process j2^J^,|^],^2|,[w| . 

Despite the generality of the problem sketched above our discussion will be focussed on fluid flow through sedi- 
mentary rocks. In particular we will discuss Fontainebleau sandstone. This model system has (together with Berea 
sandstone) acquired the status of a reference standard for modeling and analysing sedimentary rocks |p3|,p^,p|j2l]j25t| . 



General geometric characterization methods traditionally include porosities, specific surface areas, and sometimes 
correlation functions Hi ■ Recently a more refine d, qua ntitative characterization for general stochastic mi- 

crostructures was based on local porosity theory (LPT) J29|jl^ , |30| , ^l| , |l8| , ^2| -^5|] . LPT is currently the most general 
geometric characterization method because it contains as a special case also the characterization through correlation 
functions (see jl6) for details). 

Local porosity theory is used in this paper to distinguish quantitatively various models for Fontainebleau sandstone. 
More precisely, the objective of this work is to give a quantitative comparison of four microstructures. One of them 
is an experimental sample of Fontainebleau sandstone, while three of the microstructures are synthetic samples from 
computer simulation models for Fontainebleau sandstone. One of the models is a sedimentation and diagenesis model 
that tries to mimick the formation of sandstone through deposition and cementation of spherical grains. Two purely 
stochastic models generate random realizations of microstructures with prescribed porosity and correlation function. 
The first of these is based on Fourier space filtering of Gaussian random fields, and the second is based on a simulated 
annealing algorithm. 

In sectio n || we introduce and define the geometrical quantities that will be used to distinguish the microstructures. 



In section [II we present the four microstructures, their generation and characterization in terms of the generation 



procedure. In section IV we present the results and discuss the differences between the four microstructures. 



II. MEASURED QUANTITIES 



A. Porosity and Correlation Functions 



Consider a rock sample occupying a subset § C M. d of the physical space (d = 3 in the following). The sample S 
contains two disjoint subsets § = P U M with P n M = where P is the pore space and M is the rock or mineral 
matrix and is the empty set. The porosity </>(§) of such a two component porous medium is defined as the ratio 
</>(§) = V(P)/V(§) which gives the volume fraction of pore space. Here V(P) denotes the volume of the pore space, 
and V(S) is the total sample volume. 

For the sample data analysed here the set S is a rectangular parallelepiped whose sidelengths are Mi , Mi and M3 in 
units of the lattice constant a (resolution) of a simple cubic lattice. Thus the sample is represented in practice as the 
subset § = [1, Mi] x [1, M2] x [1, M3] C 1? of an infinite cubic lattice. Z denotes the set of integers, and [1, Mi] C Z 
are intervals. The position vectors = r^...^ = (fflii, ...,aid) with integers 1 < ij < Mj are used to label the lattice 
points, and r t is a shorthand notation for r ilid . A configuration (or microstructure) Z of a 2-component medium is 
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then given as 

Z=(Z 1 ,... ) Z w ) = ( Xp (r 1 ),..., Xp M) (2.1) 
where N = M1M2M3, and 

M [ 1 : for r e G , . 

*J r ) = (o : for r^G ^ 

is the characteristic (or indicator) function of a set G that indicates when a point is inside or outside of G. A stochastic 
medium is defined through the discrete probability density 

p(zi, z N ) = Prob{(Zi = zi) A ... A (Z N = z N )} (2.3) 

where Zi £ {0, 1}. Expectation values of functions f(Z) = f(z\, zjv) are defined as 
l l 

(f(zi,...,z N )) = 22-22 f{zi,...,z N )p(zi,...,z N ) (2.4) 

Zi— 3JV-0 

where the summations run over all configurations. If the medium is statistically homogeneous (stationary) then the 
average porosity is given as 



(0) = Prob{r £P} = ^ Xp (ro)) (2.5) 
where ro is an arbitrary lattice site. If the medium is also ergodic then the limit 

lim 0(S) = (0) (2.6) 

exists. There are, however, many subtleties associated with this limit (see for details). Finally, we now define the 
correlation function for a homogeneous medium as the expectation 

G(ro,r) = G(r-r )= \* m l_' m . (2.7) 
If the medium is also isotropic G(r) = G(|r|) = G(r). Obviously G(0) = 1 and G(oo) = 0. 



B. Local Porosity Distributions 

The basic idea of local porosity theory is to measure geometric observables within a bounded (compact) subset of the 
porous medium and to collect these measurements into various histograms. Let K(r, L) denote a cube of sidelength L 
centered at the lattice vector r. The set K(r, L) defines a measurement cell inside of which local geometric properties 
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such as porosity or specific internal surface are measured |30|. The local porosity in this measurement cell K(r, L) is 
defined as 



(r,L) 



V(PnK(r,£)) 
V(K(r, L) ) 



(2.8) 



where V(G) is the volume of the set G C M d . The local porosity distribution fj,(<j), L) is defined as 



M^j L) = — 5((f> — (p(r, L)) 



(2.9) 



where m is the number of placements of the measurement cell K(r, L). Ideally all measurement cells should be disjoint 
pO| ], but in practice this cannot be achieved because of poor statistics. The results presented below are obtained by 
placing K(r,L) on all lattice sites r which are at least a distance L/2 from the boundary of S, and hence in the 
following 



m 



H(Mi-L+i) 



(2.10) 



i=l 



will be used. /j,(4>,L) is the empirical probability density function (histogram) of local porosities. Its support is the 
unit interval. In the following we denote averages with respect to /i(</>, L) by an overline. Thus for a homogeneous 
and ergodic medium 



<f>(L)= / <f>ii(<j>,L)d4> = (4>) 
Jo 



(2.11) 



is the expected local porosity. In practice deviations from the last equality may occur if the measurement cells are 
overlapping. Figure |l0| below shows the average local porosity as function of L for all four samples analyzed in this 
paper showing that deviations can be as large as 0.5 percent. The deviations may be partly intrinsic and partly due 
to oversampling the central regions because the measurement cells are overlapping. Similarly the variance of local 
porosities is found as [M 



a 2 OL) = (0(L)-<ML))2 = / [4>-<KL)Yp(<l>,L)d4> 



( 



(2.12) 



where K(r , L) is any cubic measurement cell. 

It is simple to determine u(4>, L) in the limits L — > and L — > oo of small and large measurement cells. For small 
cells one finds generally |p0|.|l6| 



fi(cj>, L = 0) = ct>(&)6(<f> - 1) + (1 - <P(S)W) 



(2.13) 



where <f>(§) is the sample porosity. If the sample is macroscopically homogeneous and ergodic then one expects 



fx(4>,L ^oo) = 5(<j> - (j>(S)) 



(2.14) 
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indicating that in both limits the geometrical information contained in /i(<fi,L) consists of the single number 4>(S). 
The macroscopic limit, however, involves the questi on of macro scop ic heterogeneity versus macroscopic homogeneity 



(for more information see [|16|). In any case, if eqs. (2.13) and (2.14) hold it follows that there exists a special length 
scale L* defined as 

L* = min{£ : fi(0, L) = /i(l, L) = 0} (2.15) 
at which the (^-distributions at <p = and (f> — 1 both vanish for the first time. 

C. Local Percolation Probabilities 

The local percolation probabilities characterize the connectivity of measurement cells of a given local porosity. Let 

{1 : if K(r,L) percolates in "a" -direction 
v. • (2 ' 16) 
: otherwise 

be an indicator for percolation. What is meant by "a" -direction is summarized in Table |l]. A cell K(r, L) is called 
"percolating in the ^-direction" if there exists a path inside the set PflK(r, L) connecting those two faces of S that are 
vertical to the x-axis. Similarly for the other directions. Thus A3 = 1 indicates that the cell can be traversed along 
all 3 directions, while A c = 1 indicates that there exists at least one direction along which the block is percolating. 

The local percolation probability in the "a" -direction is now defined through 

Aa( ^ ) = ^^M (2.17) 

where <W( rj n = 1 if <fr = <?K r ! L) an d otherwise. The local percolation probability X a (<j>,L) gives the fraction of 
measurement cells of sidelength L with local porosity <f> that are percolating in the "a" -direction. 

D. Total Fraction of Percolating Cells 

The total fraction of all cells percolating along the "a" -direction is given by integration over all local porosities as 

p a (L)= f ti(<f>,L)\ a {cl>,L)d<f> (2.18) 
Jo 

This quantitiy provides an important characteristic for constructing equivalent network models. It gives the fraction 
of network elements (bond, sites etc.) which have to be permeable in an equivalent network. 

III. DESPCRIPTION OF MICROSTRUCTURES 
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A. Experimental Sample of Fontainebleau Sandstone 



The experimental sample is a threedimensional microtomographic image of Fontainebleau sandstone. This sandstone 
is a popular reference standard because of its exceptional chemical, crystallographic and microstructural simplicity 
^3j,p4|. Fontainebleau sandstone consists of monocrystalline quartz grains that have been eroded for long periods 
before being deposited in dunes along the shore during the Oligocene, i.e. roughly 30 million years ago. It is well sorted 
containing grains of around 200/xm in diameter. During its geological evolution, that is still not fully understood, the 
sand was cemented by silica crytallizing around the grains. Fontainebleau sandstone exhibits intergranular porosity 
ranging from 0.03 to roughly 0.3 p4| . 

The computer assisted microtomography was carried out on a micro-plug drilled from a larger original core. The 
original core from which the micro-plug was taken had a porosity of 0.1484, a permability of 1.3-D and a formation 
factor of 22.1. The porosity </>(§ex) of our microtomographic data set is only 0.1355(see Table ||). The difference 
between the porosity of the original core and that of the final data set is due to the heterogeneity of the sandstone 
and to the difference in sample size. The experimental sample is referred to as EX in the following. The pore space 
of the experimental sample is visualized in Figure [l]. 



B. Sedimentation, Compaction and Diagenesis Model 



The sedimentation, compaction and diagenesis model, abbreviated as DM in the following, is obtained by 
numerically modelling the main geological sandstone- forming processes 28 ll|. Image analysis of backscattered 



electron/cathodo-lumincsccnce images of thin sections provides input data such as porosity, grain size distribution, 
a visual estimate of the degree of compaction, the amount of quartz cement and clay contents and texture. The 
sandstone modelling is carried out in three main steps: grain sedimentation, compaction and diagenesis. Here we give 



only a rough sketch of the algorithms and refer the reader to 28 15J for a detailed description 



Grain sedimentation commences with image analysis of thin sections. The grain size distribution is measured using 
an erosion-dilation algorithm. Spherical grains with random diameters chosen from the grain size distribution are 
dropped onto the grain bed and relaxed into a potential energy minimum. The sedimentation environment may be 
low-energy (local minimum) or high-energy (global minimum). 

Compaction reduces the bulk volume (and porosity) in response to vertical stress from the overburden. It is 
modelled here as a linear process in which the vertical coordinate of every sandgrain is shifted vertically downwards 
by an amount proportional to the original vertical position. The proportionality constant is called the compaction 
factor. Its value for our Fontainebleau sandstone is estimated to be 0.1 from thin section analysis. 

In the diagenesis part only a subset of known diagenetical processes are simulated, namely quartz cement overgrowth 
and precipitation of authigenic clay on the free surface. Quartz cement overgrowth is modeled by radially enlarging 
each grain. If Rq denotes the radius of the originally deposited spherical grain, its new radius along the direction r 
from grain center is taken to be [7|Jl5t| 

R(r) = R + mm(a£{r)i , £{r)) (3.1) 

where t(r) is the distance between the surface of the original spherical grain and the surface of its Voronoi ployhedron 
along the direction r. The constant a controls the amount of cement, and the growth exponent 7 controls the type of 
cement overgrowth. For 7 > the cement grows preferentially into the pore bodies, for 7 = it grows concentrically, 
and for 7 < quartz cement grows towards the pore throats |15|| . Authigenic clay growth is simulated by precipitating 
clay voxels on the free mineral surface. The clay texture may be pore-lining or pore-filling or a combination of the 
two. 

For modeling the Fontainebleau sandstone we used a compaction factor of 0.1, and the cementation parameters 
7 = —0.6 and a = 2.9157. The resulting configuration of our sample DM is displayed in Figure g. 
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C. Gaussian Field Reconstruction Model 



The Gaussian field (GF) reconstruction model provides a random pore space configuration in such a way that its 
correlation function Ggf(i") equals a prescribed reference correlation function Gq(t). In our case Go(r) = Gex(i") 
the reference is the correlation function of the experimental sample described above. The method of Gaussian field 
reconstruction is well documented in the literature |2(]|j^|]|37| , and we shall only make a few remarks that the reader 
may find of interest when implementing the method. 

Given the reference correlation function Gex(i") and porosity 0(§ex) of the experimental sample the three main 
steps of constructing the sample §gf with correlation function Gqf(t) = Gex(i") are as follows: 

1. A standard Gaussian field X(r) is generated which consists of statistically independent Gaussian random vari- 
ables Xgtat each lattice point r. 

2. The field X(r) is first passed through a linear filter which produces a correlated Gausssian field Y(r) with 
zero mean and unit variance. The reference correlation function Gex(i") and porosity </>(§ex) enter into the 
mathematical construction of this linear filter. 

3. The correlated field Y (r) is then passed through a nonlinear discretization filter which produces the reconstructed 
sample §gf- 

Details of these three main steps are documented in Ref. [^,^6|. However, in these traditional methods the process 
described in step || is computationally difficult because it requires the solution of a very large set of non-linear 
equations. We have followed an alternate and computationally more efficient method proposed in Ref. S which 
uses Fourier Transforms. For the sake of completeness we briefly describe this. Later we shall discuss some of the 
difficulties experienced while implementing this. 

In the Fourier transform method the linear filter in step || is defined in Fourier space through 

F(k)=a(G y (k))ix(k), (3.2) 
where M = M\ = M% = M3 is the sidelength of a cubic sample, a = Mi is the normalisation factor, and 

r 

denotes the Fourier transform of X(r). Similarly Y(k) is the Fourier transform of Y(r), and Gy(k) is the Fourier 
transform of the correlation function Gy(r). Gy(r) has to be computed by an inverse process from the correlation 
function Gex(i") and porosity of the experimental reference (details in ||). 

It is important to note that the Gaussian field reconstruction requires a large separation £ex -C iV 1 ^ where £ex is 
the correlation length of the experimental reference, and N = M1M2-M3 is the number of sites. £ex is defined as the 
length such that Gex(^) ~ for r > £ex- If the condition £ex -C N 1 ^ is violated then step || of the reconstruction 
fails in the sense that the correlated Gaussian field Y(r) does not have zero mean and unit variance. In such a 
situation the filter Gy(k) will differ from the Fourier transform of the correlation function of the Y(r). It is also 
difficult to calculate Gy(r) accurately near r = |. This leads to a discrepancy at small r between Ggf(^) and 
GexM- The problem can be overcome by choosing large M as we verified in d = 1 and d — 2. However, in d = 3 
very large M also demands prohibitively large memory. In earlier work (36]|| the correlation function Gex(i") was 
sampled down to a lower resolution, and the reconstruction algorithm then proceeded with such a rescaled correlation 
function. This leads to a reconstructed sample §gf which also has a lower resolution. Such reconstructions have 
lower average connectivity compared to the original model |38| Because we intend a quantitative comparison with the 
microstructure of Sex it is necessary to retain the same level of resolution. Hence we use throughout this article the 
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original correlation function Gex(t) without subsampling. Because GexM is nearly for r > 30a we have truncated 
Gex(^) at r = 30a to save computer time. The final configuration Sgf with M = 256 generated by Gaussian filtering 
reconstruction that is used in the comparison to experiment is displayed in Figure 0. 



D. Simulated Annealing Reconstruction Model 

The simulated annealing (SA) reconstruction model is a second method to generate a threedimcnsional random 
microstructurc with prescribed porosity and correlation function. A simplified implementation was recently discussed 
m Ref . [|IJ and we follow their algorithm here. The method generates a configuration §sa by minimizing the deviations 
between Gsa(i") and a predefined reference function Gq(t). Of course in our case we have again the Fontainebleau 
sandstone as reference, i.e. Go(r) = Gex(f). 

The reconstruction is performed on a cubic lattice with side length M = Mi = M2 = M3 and lattice spacing a. 
The lattice is initialized randomly with 0's and l's such that the volume fraction of 0's equals </>(§ex)- This porosity 
is preserved throughout the simulation. For the sake of numerical efficiency the autocorrelation function is evaluated 
in a simplified form using ]2lf| 

G SA (r) (Gsa(0) - Gsa(O) 2 ) +G sa (0) 2 = 

= * M ( r ) (x M ( r + rei ) + V r + re2 ) + * M ( r + re3 )) ( 3 - 4 ) 

r 

where e.; are the unit vectors in direction of the coordinate axes, r = 0, . . . , M- — 1, and where a tilde ~ is used to 
indicate the directional restriction. The sum ^ r runs over all M 3 lattice sites r with periodic boundary conditions, 
i.e. Ti + r is evaluated modulo M. 

We now perform a simulated annealing algorithm to minimize the "energy" function 

£ = ^(GsA(r)-G E x(r)) 2 , (3.5) 



defined as the sum of the squared deviations of G$a from the experimental correlation function Gex- Each update 
starts with the exchange of two pixels, one from pore space, one from matrix space. Let n denote the number of 
the proposed update step. Introducing an acceptance parameter T„, which may be interpreted as an n-dependent 
temperature, the proposed configuration is accepted with probability 

v = min y 1 , ex p y- Er T E En ^ 1 J ^ • (3-6) 

Here the energy and the correlation function of the configuration is denoted as E n and GsA,n, respectively. The 
evaluation of Gsa.m does not require a complete recalculation. It suffices to u pdate the correlation function GsA,n-i 
of the previous configuration by adding or subtracting those products in (^J) which changed due to the exchange of 
pixels. In case the proposed move is rejected, the old configuration is restored. 

The generation of a configuration with correlation Gex is achieved by lowering T. At low T the system approaches 
a configuration that minimizes the energy function. In our simulations we lower T n with n as 

T « = exp (-IobW>- (3J) 
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We stop the simulation when 20000 consecutive updates are rejected. In our simulation this happened after 2.5 x 10 8 
updates (~ 15 steps per site). The resulting configuration §sa for the simulated annealing reconstruction is displayed 
in Figure ||. 



Our definition of the correl atio n function in (3.4) deserves some comment. A complete evaluation of the corre- 
lation function as defined in (2.7) requires such a great numerical expense that the algorithm is too slow to allow 
threcdimensional reconstructions within a reasonable time. Therefore, to increase the speed of the algorithm, the 
correlation function is only evaluated along the directions of the coordinate axes as indicated in (3.4). As a result of 
this simplification the reconstructed sample may cease to be isotropic. It will in general deviate from the reference 
correlation function in all directions other than those of the axes. In the special case of the correlation function of 
the Fontainebleau san dsto ne, however, this effect seems to be small (see below). This may serve as an a posteriori 
justification for using (3.4). 



IV. RESULTS AND DISCUSSION 



We begin our presentation of the results with an analysis of traditional quantities such as porosities and correlation 
functions of the four samples. Then we proceed to a visual characterization of the threedimcnsional images. Next 
we shall discuss local porosities and percolation probabilities, and finally we conclude with implications for transport 
properties. 



A. Conventional Analysis 



Table | gives a synopsis of different properties of the four samples. The preparation of the various samples was 
described in detail in section [II. The dimensions and porosities also need no further comment. Samples GF and 
SA were constructed to have the same correlation function as sample EX. This is indicated in the line labeled G(r). 
In Figure || we plot the directionally averaged correlation functions G(r) = (G(r, 0, 0) + G(0, r, 0) + G(0, 0, r))/3 of the 
four samples where G(ri,r2,rs) = G(r). Gdm(^) differs clearly from the rest. Acci dentally, however, Gdm(0, 0,r) w 
Gex(0, 0,r). Gcf(r) differs from Gex(^) for small r as discussed in section III C| above. Remember also that by 



construction GgfM is not expected to equal Gex(i") for r larger than 30. The discrepancy at small r reflects the 
quality of the linear filter, and it is also responsible for the differences of the porosity and specific internal surface. 
Although the reconstruction method of sample §sa is intrinsically anisotropic the correlation function of sample 
SA agrees also in the diagonal directions with that of sample EX. Sample §dm on the other hand has an anisotropic 
correlation function. 

If two samples have the same correlation function they are also expected to have also the same specific internal 
surface as calculated from 



S=-4(0)(l-(0)) dG(r) 



dr 



(4.1) 



The next line in Table || labeled S gives the specific internal surfaces. 

If one defines a decay length by the first zero of the correlation function then the decay length is roughly 18a for 
samples EX, GF and SA. For sample DM it is somewhat smaller mainly in the x- and y-direction. The correlation 
length, which will be of the order of the decay length, is thus relatively large compared to the system size. Together 
with the fact that the percolation threshold for continuum systems is typically around 0.15 this might explain why 
models GF and SA are connected in spite of their low value of the porosity. 

In summary, the samples §gf and Ssa were constructed to be indistinguishable with respect to porosity and 
correlations from §ex- The imperfection of the reconstruction method for sample GF , however, accounts for the 
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deviations of its correlation function at small r from that of sample EX. 



B. Visual inspection of images 



We now collect results from a visual comparison. Visual inspection of Figures |lj through ^ reveals that none of the 
models §dm,§gf or §sa resemble closely the experimental microstructure §ex- This applies in particular to samples 
GF and SA which were constructed to match the traditional geometrical characteristics of sample EX, such as porosity, 
specific surface and correlation function. 

Figures [j] through || suggest that samples §gf and Ssa have isolated islands of matrix space although this cannot be 
seen directly because the pore space is rendered opaque. Isolated islands of matrix space cannot exist in a real porous 
medium such as sample EX. They are also absent in the compaction and diagenesis model DM. The comparison is 
indicated in the line labeled "isolated M" in Table The pore surfaces in samples GF and SA are much rougher 
than in samples EX and DM. Sample DM appears visually more homogeneous than the other samples. Although 
there is no anisotropy visible for sample DM from Figure |^ its connectivity properties will be found below to be 
strongly anisotropic. 

In summary the traditional characteristics such as porosity, specific surface and correlation functions are insufficient 
to distinguish different microstructures. Visual inspection of the pore space by the human eye indicates that samples 
GF and SA have a similar structure which, however, differs from the structure of sample EX. Although sample 
DM resembles sample EX more closely with respect to surface roughness it differs visibly in the shape of the grains. 



C. Local Porosity Analysis 



Wo turn to an analysis of the fluctuations in local porosities. The differences in visual appearance of the microstruc- 
tures find a quantitative expression here. 

The local porosity distributions fi((j), 20) of the four samples at L — 20a are displayed as the solid lines in Figures 
H through ||. The ordinates for these curves are plotted on the right vertical axis. The figures show that the original 
sample exhibits stronger porosity fluctuations than the three model samples except for sample SA which comes close. 
Sample DM has the narrowest distribution which indicates that it is most homogeneous. Figures ^-|| show also that 
the component at the origin, //(0,20), is largest for sample EX, and smallest for sample GF. For samples DM and 
SA the values of /i(0, 20) are intermediate and comparable. Plotting /i(0, L) as a function of L we find that this 
remains true for all L. These results indicate that the experimental sample EX is more strongly heterogeneous than 
the models, and that large regions of matrix space occur more frequently in sample EX. A similar conclusion may be 
drawn from the variance of local porosity fluctuations which will be studied below. The conclusion is also consistent 
with the results for L* shown in Table |2|. L* gives the sidelength of the largest cube that can be fit into matrix space, 
and thus L* may be viewed as a measure for the size of the largest grain. Table ^ shows that the experimental sample 
has a larger L* than all the models. It is interesting to note that plotting (j,(l,L) versus L also shows that the curve 
for the experimental sample lies above those for the other samples for all L. Thus, also the size of the largest pore 
and the pore space heterogeneity are largest for sample EX. If /j,(cf),L*) is plotted for all four samples one finds two 
groups. The first group is formed by samples EX and DM, the second by samples GF and SA. Within each group 
the curves /ii(</>,L*) nearly overlap, but they differ strongly between them. 

Figures |l(], [ll], and |l2| exhibit the dependence of the local porosity fluctuations on L. In Figure |ll] we plot the 



variance of the local porosity fluctuations, defined in eq.(2.12) as function of L. The variances for all samples indicate 
an approach to a (5-distribution according to eq. ( 2.14| ). Again sample DM is most homogeneous in the sense that 



its variance is smallest. The agreement between samples EX, GF and SA reflects the agreement of their correlation 







functions, and is expected by virtue of eq. (2.12). Figure [12] shows the skewness as a function of L calculated from 



K3{L) = — — (4 - 2) 



where cr(L) is the variance defined in eq. (2.12). K3 characterizes the asymmetry of the distribution, and the difference 
between the most probable local porosity and its average. Again samples GF and SA behave similarly, but sample 
DM and sample EX differ from each other, and from the rest. 

At L = 4a the local porosity distributions n(4>, 4) show small spikes at equidistantly spaced porosities for samples 
EX and DM, but not for samples GF and SA. The spikes indicate that models EX and DM have a smoother surface 
than models GF and SA. For smooth surfaces and small measurement cell size porosities corresponding to an interface 
intersecting the measurement cell occur with higher frequency, and this gives rise to spikes. The presence of isolated 
islands of pore or matrix space reduces these spikes. It is unclear at present whether the spikes persist when the 
measurement cells are chosen to be nonoverlapping. 



D. Local Percolation Analysis 



Visual inspection of Figures [j] through |] did not allow us to recognize the degree of connectivity of the various 
samples. A quantitative characterization of the connectivity is provided by the local percolation probabilities p0| , ^5[ , 
and it is here that the samples differ most dramatically. 

The samples EX, DM , GF and SA are globally connected in all three directions. This, however, does not imply 
that they have similar connectivity. The last line in Table || gives the fraction of blocking cells at the porosity 0.1355 
and for L*. It gives a first indication that the connectivity of samples DM and GF is, in fact, much poorer than that 
of the experimental sample EX. 

Figures ^ through [9] give a more complete account of the situation by exhibiting X a (</>, 20) for a — 3, c, x, y, z for all 
four samples. First one notes that sample DM is strongly anisotropic in its connectivity. It has a higher connectivity 
in the z-direction than in the x- or j/-direction. This might be due to the anisotropic compaction process. X z (<p, 20) 
for sample DM differs from that of sample EX although their correlation functions in the z-direction are very similar. 
The A-functions for samples EX and DM rise much more rapidly than those for samples GF and SA. The inflection 
point of the A-curves for samples EX and DM is much closer to the most probable porosity (peak) than in samples 
GF and SA. All of this indicates that connectivity in cells with low porosity is higher for samples EX and DM than 
for samples GF and SA. In samples GF and SA only cells with high porosity are percolating on average. In sample 
DM the curves A^, Aj, and A3 show strong fluctuations for A « 1 at values of <fi much larger than the (</)) or 0(Sdm)- 
This indicates a large number of high porosity cells which are nevertheless blocked. The reason for this is perhaps 
that the linear compaction process in the underlying model blocks horizontal pore throats and decreases horizontal 
spatial continuity more effectively than in the vertical direction, as shown in p8[ ], Table 1 p. 142. 

The absence of spikes in fi(<f>,4) for samples GF and SA combined with the fact that cells with average porosity 
(« 0.135) are rarely percolating suggests that these samples have a random morphology similar to percolation. 



E. Implications for Transport Properties 



The connectivity analysis of local porosity theory allows to make some predictions for transport properties (such 
as conductivity or permeability) without actually calculating them. A detailed comparison between the predictions 
of local porsity theory and exact calculation of transport proper ties w ill appear elsewhere [ p9| These predictions are 
made by calculating the total fraction of percolating cells eq. (2.18). The insets in Figures rc] through show the 
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functions p a (L) = X a (<p,L) for a = 3,%,y,z,c for each sample. The curves for samples EX and DM are similar 
but differ from those for samples GF and SA. In Figure |l3| we plot the curves p^{L) of all four samples in a single 
figure. The samples fall into two groups {EX, DM} and {GF,SA} that behave very differently. Figure [l^ shows that 
reconstruction methods based on correlation functions do not reproduce the connectivity properties of porous 

media. As a consequence, within the effective medium approximation of local porosity theory [ |30[ samples GF and 
SA would both yield much lower permeabilities or conductivities than those of samples EX and DM. Based on 
these results it appears questionable whether correlation function reconstruction can produce reliable models for the 
prediction of transport. 
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TABLE CAPTIONS 



Table 1: Legend for index a of local percolation probabilities \ a {<j>, L). 
Table 2: Overview of various properties for the four samples 
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TABLE § 



index a 


meaning 


X 


ir-direction 


y 


y-direction 


z 


z-direction 


3 


(x A y A z)-direction 


c 


(x V y V z)-direction 



TABLE | 



Properties 


§ex 


Sdm 




S S A 


Origin 


Experiment 


Diagenesis Model 


Gaussian Field 


Simulated Annealing 


Mx x M 2 x M 3 


300 x 300 x 299 


255 x 255 x 255 


256 x 256 x 256 


256 x 256 x 256 


0(S) 


0.1355 


0.1356 


0.1421 


0.1354 


G(r) 


Gex 


Gdm 


Ggf ~ Gex 


Gsa = Gex 


^from f | r=0 


0.078 


0.082 


0.125 


0.083 


Isotropy 


xyz 


xy 


xyz 


xyz 


isolated M 


No 


No 


Yes 


Yes 


Pore surface 


smooth 


smooth 


rough 


rough 


L* 


35 


25 


23 


27 


Connectivity 


xyz 


xyz 


xyz 


z 


1 - A C (0.1355,L*) 


0.0045 


0.0239 


0.3368 


0.3527 
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FIGURE CAPTIONS 



Figure 1: Threedimensional pore space of Fontainebleau sandstone (sample EX). The resolution of the image 
is a — 7.5/xm, the sample dimensions are Mi — 300, M 2 = 300, M 3 = 299. The porosity is 0(§ex) = 
0.1355. The pore space is indicated opaque, the matrix space is transparent. The lower image shows 
the front plane of the sample as a twodimcnsional thin section (pore space black, matrix grey). 

Figure 2: Threedimensional pore space of the sedimentation and diagenesis model (sample DM). The resolution 
is a = 7.5/zm, the sample dimensions are Mi — 255, M 2 = 255, M 3 = 255. The bulk porosity is 
0(§dm) = 0.1356. The pore space is indicated opaque, the matrix space is transparent. The lower 
image shows the front plane of the sample as a twodimensional thin section (pore space black, matrix 
grey). 

Figure 3: Threedimensional pore space having the same correlation function as the experimental sample of 
Fontainebleau sandstone (sample GF). The pore space was constructed using Gaussian random fields 
which are subsequently filtered. The resolution is a = 7.5/im, the sample dimensions are Mi = 256, 
Mi = 256, M3 = 256. The bulk porosity is </>(§>gf) = 0.1421. The pore space is indicated opaque, the 
matrix space is transparent. The lower image shows the front plane of the sample as a twodimensional 
thin section (pore space black, matrix grey). 

Figure 4: Threedimensional pore space having the same correlation function as the experimental sample of 
Fontainebleau sandstone (sample SA). The pore space was constructed using a simulated annealing 
algorithm. The resolution is a = 7.5/zm, the sample dimensions are Mi = 256, M 2 = 256, M 3 = 256. 
The bulk porosity is </>(§sa) = 0.1354. The pore space is indicated opaque, the matrix space is 
transparent. The lower image shows the front plane of the sample as a twodimensional thin section 
(pore space black, matrix grey). 

Figure 5: Averaged directional correlation functions G(r) = (G(r, 0,0) + G(0,r, 0) + G(0,0, r))/3 of all four 
samples. 

Figure 6: Local percolation probabilities A Q (</>, 20) (broken curves, values on left axis) and local porosity dis- 
tribution 20) (solid curve, values on right axis) at L = 20 for sample EX. The inset shows the 
function p a (L). The line styles correponding to a = c, x, y, z, 3 are indicated in the legend. 

Figure 7: Local percolation probabilities \ a {<p, 20) (broken curves, values on left axis) and local porosity dis- 
tribution /j,(cj>, 20) (solid curve, values on right axis) at L — 20 for sample DM. The inset shows the 
function p a (L). The line styles correponding to a = c, x, y, z, 3 are indicated in the legend. 

Figure 8: Local percolation probabilities X a (<p, 20) (broken curves, values on left axis) and local porosity dis- 
tribution n((j>, 20) (solid curve, values on right axis) at L = 20 for sample GF. The inset shows the 
function p a (L). The line styles correponding to a = c, x, y, z, 3 are indicated in the legend. 

Figure 9: Local percolation probabilities \ a (cf), 20) (broken curves, values on left axis) and local porosity dis- 
tribution fi((j), 20) (solid curve, values on right axis) at L — 20 for sample SA. The inset shows the 
function p a (L). The line styles correponding to a = c,x,y, z,3 are indicated in the legend. 

Figure 10: Average local porosities for sample EX(solid line with tick) DM (dashed line with cross) GF(dotted 
line with square), and SA(dash-dotted line with circle). 

Figure 11: Variance of local porosities for sample EX (solid line with tick) DM (dashed line with cross) GF (dotted 
line with square), and SA(dash-dotted line with circle). 

Figure 12: Skewness of local porosities for sample EX(solid line with tick) DM (dashed line with cross) GF(dotted 

line with square), and SA(dash-dotted line with circle). 
Figure 13: p^{L) for sample EX(solid line with tick) DM (dashed line with cross) GF(dotted line with square), 

and SA(dash-dotted line with circle). 
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Figure |ll| 
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Figure |l3| 
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